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SUMMARY 


An algorithm for maximum likelihood (ML) estimation is developed with an effi- 
cient method for approximating the sensitivities. The algorithm is applicable to any 
parameter -estimation problem and is particularly suited for multivariable dynamic 
systems. The ML algorithm relies on a new optimization method closely related to a 
modified Newton -Raphs on (MNR) technique; the new optimization method is referred to 
as a modified Newton-Raphson with estimated sensitivities (MNRES ) • 

MNRES determines sensitivities by using slope information from local surface 
approximations of each output variable in parameter space. The fitted surface allows 
sensitivity information to be updated at each iteration with a significant reduction 
in computational effort. MNRES determines the sensitivities with less computational 
effort than using either a finite-difference method or integrating the analytically 
determined sensitivity equations as in an MNR procedure. The choice of the type of 
surface (for example, nth -order polynomial or spline) and the method of fitting the 
surface (for example, least squares or simply solving simultaneous equations) is made 
by the user to suit the particular need. 

Two surface-fitting methods are discussed and demonstrated, while other possi- 
bilities are indicated. Comparisons are made between MNRES and other commonly used 
optimization methods such as the Nelder— Mead simplex method (a search technique) and 
the modified Newton-Raphson method (a gradient technique). Several sample problems 
are solved to compare the techniques. Simple linear systems are used at first, and 
then nonlinear aircraft estimation problems are solved by using both real and simu- 
lated data. MNRES is found to be equally accurate and substantially faster than the 
commonly used techniques. The reduction in computational effort provided by MNRES is 
dependent on several factors: the choice of surface-fitting method, the number of 

unknown parameters, data quality, and accuracy of the sensitivity calculations. 

MNRES eliminates the need to derive sensitivity equations for each model, thus pro- 
viding flexibility to use model equations in any format that is convenient and pro- 
viding a more generally applicable algorithm. 


INTRODUCTION 

System-identification and parameter-estimation techniques used to determine 
aircraft stability and control parameters from flight data are well established for 
the linear flight regimes. In these flight regimes, the aircraft aerodynamic model 
can be expressed as a linear function of the states and control inputs. In the non- 
linear regimes, as in high-angle-of-attack conditions, the aerodynamic model can be 
a complex function of states and control inputs. Thus, identifying the best mathe- 
representation (model structure) and estimating model parameters can be very 
difficult. The need to improve identification and estimation techniques for non- 
linear aircraft models has led to further development of the methodology for system 
identification. 


In reference 1 a modified stepwise regression (MSR), along with several testing 
criteria, is suggested as an efficient method to determine model structure and obtain 
initial parameter estimates. This was found to be a useful technique for high-angle- 
of-attack aircraft-identification problems in which nonlinear aerodynamic effects are 



present. However, the estimates are asymptotically biased and their standard errors 
are based on a simplified expression which is valid only for the "classical" linear 
regression. Therefore, it is beneficial to improve the MSR estimates by using the 
maximum likelihood (ML) method which has more favorable asymptotic properties 
(ref. 2). In addition, it is desirable for the ML algorithm to be independent of 
model structure and sensitivity calculations; this independence provides a more 
generally applicable algorithm and potentially eliminates some of the computational 
burden found with integrating sensitivity equations. Also, a very efficient ML esti- 
mation algorithm is desirable to reduce the computational effort involved in process- 
ing a large number of parameters and candidate models. Reducing computational 
requirements of the ML method requires careful examination of the optimization 
methods utilized in the algorithm. Although nonlinear, unconstrained optimization 
problems have been studied quite extensively (ref. 3), little has been done to 
improve the optimization techniques as they apply to aircraft estimation problems. 

The numerical aspects of computing ML estimates for linear dynamic systems in state- 
vector form was studied in reference 4. 

Aircraft estimation problems, which belong to a class of problems involving 
dynamic systems, require substantial computational effort at each step of the opti- 
mization process. At each step the parameters are updated and the equations of 
motion are integrated to obtain time histories of each response variable. Many ML 
estimation procedures for aircraft problems use a Newton-Raphson optimization method 
(ref. 5) that requires solving (integrating) sensitivity equations. This accounts 
for most of the computational effort since an m state system with n unknown 
parameters requires mn sensitivity equations plus m state equations to be inte- 
grated each iteration. Several states and 30 or 40 parameters are not unrealistic 
for aircraft estimation problems. This is for only one flight condition; if a model 
is desired throughout the entire flight envelope, the computational requirements 
become overwhelming since analysis of various flight conditions may require more than 
one candidate model. 

The objective of this report is to provide an ML method which does not require 
the analytical form of the sensitivity equations to formulate the algorithm and which 
provides a method with reduced computational requirements compared with the commonly 
used algorithms. The report first considers both gradient and search optimization 
methods: the Nelder-Mead simplex method (SM) (a search technique (ref. 6)) and the 

commonly used modified Newton-Raphson method (MNR) (a gradient technique). After a 
comparison of these methods, the MNR is selected as a method for further considera- 
tion. As a main objective of the report, an efficient algorithm is subsequently 
proposed. This algorithm will be referred to as a modified Newton-Raphson with 
estimated sensitivities (MNRES ) method. It determines sensitivities without the 
analytical form of the sensitivity equations and does it more efficiently than a 
finite-difference method. Therefore, it is easily used with any model structure 
determined, for example, by the MSR method. Properties of the improved method are 
discussed and several test cases are provided. Finally, the MNRES algorithm is 
applied to a nonlinear aircraft estimation problem by using both real and simulated 
data. For this study the estimation problem is simplified by assuming that no 
process noise exists in the system dynamics. 
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system matrix for state equation 
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control-distribution matrix for state equation 
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rolling-moment coefficient, M /qSb 
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number of elements in B matrix 

least-squares -parameter covariance matrix divided by a 
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dynamic pressure, — pv , Pa 
measurement-noise covariance matrix 
yaw rate, rad/sec 
wing area, m 2 

vector of sensitivities for kth element of Y 
sensitivity of kth element of Y to Xth element of 0 
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processing time for Control Data Corporation CYBER 175 digital 
computer, sec 

input vector 

airspeed, m/sec 

vector of measurement noise at ith data point 
vector of states or matrix of parameters in MNRES 
vector of outputs 
scalar output 
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incremental value 
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Kronecker delta 
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v. vector of residuals at ith data point 

1 

3 

p air density, kg/m 

a standard error 

<|> roll angle, rad 

Subscripts: 

a partial derivative with respect to indicated quantity 

H,L highest and lowest cost, respectively 

i,j,k,Jl general indices 
Superscript: 

j index of surface-fitting points 

Abbreviations: 

CPU central processing unit 

MAX ML program using MNRES algorithm 

MAXLIK ML program using MNR algorithm 

ML maximum likelihood 

MNR modified Newton -Raphs on 

MNRES modified Newton -Raphs on with estimated sensitivities 

MSR modified stepwise regression 

NR Newton-Raphson 

SM simplex method 

Matrix exponents: 

T transpose matrix 

-1 inverse matrix 

Mathematical notation: 

* estimated quantity when over symbol 

. derivative with respect to time when over symbol 

V gradient operator 
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COMPARISON OF SEARCH AND GRADIENT OPTIMIZATION METHODS 


The twofold objective of this report is to develop an optimization method for 
use in the ML estimation which (1) eliminates the need to solve sensitivity equations 
and (2) improves, if possible, the speed of the overall estimation algorithm. This 
should be accomplished while considering practical examples of nonlinear, multi- 
variable dynamic systems such as aircraft systems, which are of primary interest for 
this research. As a result of the twofold objective for this report, two well- 
established optimization methods were chosen for consideration. 

Two optimization methods (the Nelder-Mead simplex method (SM), also called the 
flexible polyhedron method, and the modified Newton-Raphson (MNR) method) will be 
described in this section. More details are provided in references 3 and 6 for the 
SM and in references 2 and 5 for the MNR. A variation of the MNR will be used in 
which the derivative information (sensitivities) is computed by using finite differ- 
ences (refs. 7 and 8). Both forms of the MNR are used in this report. 

The Nelder-Mead algorithm was chosen for this study because search methods ini- 
tially appear to be good candidates for reducing computational demands in aircraft 
estimation problems. These methods avoid derivative calculations, which is where 
most of the computation time is spent in the quasi-Newton methods. The search meth- 
ods are independent of model form and thus are readily applicable to any aerodynamic 
model. The SM is used at the Langley Research Center and has been found to be advan- 
tageous in some aircraft design and control applications (ref. 9). 

This report is primarily concerned with nonlinear aircraft estimation problems. 
Since the MNR approach is commonly used for these problems, it is included as a 
benchmark. Although it is computationally burdensome to estimate derivatives, this 
information enables relatively fast convergence of the optimization methods. In 
fact, Newton's method converges in one pass for cost functions which are quadratic. 
Hence, Newton-Raphson gradient techniques used for estimation problems of dynamic 
systems are expected to converge faster when the quadratic approximations for the 
cost functions are good. Also Newton's method and the quasi-Newton methods provide 
both step size and direction for each iteration independent of relative parameter 
scaling. In some problems, however, additional control of step size is needed to 
ensure convergence. Since removing the requirement to solve sensitivity equations 
has been proposed, the MNR algorithm in this report will use a simple finite- 
difference method .except when otherwise noted. This does not turn out to be too 
costly in terms of computational time (refs. 7 and 8). However, care must be taken 
to obtain the derivatives as accurately as possible, especially in large problems. 


Description of Two Candidate Optimization Methods 

The SM represents a class of methods, called search methods, which do not use 
derivative information to minimize a scalar function of n variables J(0). The 
method uses a flexible polyhedron surface with n + 1 vertices where each vertex is 
defined by a vector 0. The vertex 0 H , producing the highest value of J(0), is 
projected through the centroid of the remaining vertices to define a new vertex. 

This new vertex, and the remaining ones without 0„, form a new polyhedron. This 
operation is called a "reflection." If the new vertex produces a lower cost than 
0 L (the vertex producing the smallest J(0)), then an expansion takes place and a 
new vertex is located farther out along the same projection. Similarly, if higher 
costs are found, a contraction takes place. The minimum of the cost function is 
found by repeatedly deleting the point having the highest value of J(0) and adding 
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new projected points that produce lower J(9). The flexible polyhedron is able to 
adapt to the shape of J(0) by stretching down slopes, contracting near minima, and 
changing direction in curved valleys. 

The modified Newton-Raphson optimization method belongs to a class of methods 
known as variable— metric or quasi— Newton methods. This class has various techniques 
for approximating the Hessian matrix of the cost function (matrix of second partials 
with respect to parameters), but they all use first-derivative information for the 
approximation. The MNR method will be illustrated in an ML estimation problem for a 
linear system without process noise. The theory behind ML will not be discussed; see 
reference 2 for details. It is assumed that only the measured outputs are corrupted 
by noise and that the noise is an uncorrelated, Gaussian white-noise source with zero 
mean. The problem is to minimize the errors between the computed model outputs and 
the actual measured outputs. The optimization problem is nonlinear and MNR is 
usually applied because of its good convergence characteristics even for large num- 
bers of unknown parameters. The system equations are as follows: 

State equation: 

X = AX + BU X(0) = X Q (1 ) 

Output equation: 

Y = CX + DU (2) 

Measurement equation: 

Z i = Y i + V i (i = 1, 2, ..., N) (3) 

where it is assumed that 



The symbol R represents a diagonal measurement-noise covariance matrix. Without 
process noise, the cost function to be minimized is 


J(0) 


- 1 £ (h - ? i) v ’( z i - ? i) * I ln l R ''l 


i=1 


(5) 


The unknown R can be determined by minimization of the cost with respect to R. 
This minimization produces (see ref. 2) 


1 ^ T 

R -« £ VV 

1=1 


(6) 


8 



The cost can then be written as 


^ i.1 A A 

= T Z ( Z i ~ ) T k 1 (z. - Y. ) + Constant ( 7 ) 

i=1 11 


This is the same cost function used in the output error technique from reference 2, 
except that the covariance of the measurement noise is used as a weighting matrix. 

The MNR method accomplishes the minimization by expanding Y (the computed output 
vector) about 0 Q (the initial unknown parameter vector). A Taylor series expansion 
truncated to first order is computed as 


A9 (8) 


where A0 =0 - 0 Q . Then, by substituting equation (8) into (7), a quadratic 
approximation of J(0) is obtained. The increment A0 is the unknown. Differenti- 
ating J with respect to 0 and equating the derivative to 0 to find the minima 
results in 
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This is often written as 
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(9) 


( 10 ) 


(11 ) 


( 12 ) 


emphasizing the Fisher information matrix M and gradient terms. So, for each kth 
iteration the new estimate 0 fc+1 is given as 0 k+i = 0 k + A0 k+1 . Convergence is 
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achieved when AJ/J and A9/9 are small enough. The sensitivities are deter- 

mined separately from the aforementioned steps. This can be done numerically with a 
simple finite-difference method or by integrating the sensitivity equations. The 
sensitivity equations for the system defined in equations (1) and (2) are 


X = AX + X 

a- a. oa. 

3 3 3 


(j “ 1 f * . * f n_ ) 


(13) 


• 5B 

x b. = ^b. + ab; 
3 l 1 


u 


(j — If 2f ...f n b ) 


(14) 


where each element of A and B is an unknown parameter. Details of the 
sensitivity-equation approach are discussed in reference 2 and the finite-difference 
approach is discussed in references 7 and 8. 


Performance of Methods on Estimation Problem I 

Estimation problem I demonstrates and compares the Nelder-Mead and modified 
Newton-Raphson optimization methods in a simple ML estimation problem. The MNR 
method uses a finite-difference technique to compute derivatives which satisfies the 
first objective of this research, eliminating the need to compute analytical gra- 
dients. The MNR method using numerically determined derivatives generally performs 
with about the same speed as a method using analytically determined derivatives, that 
is, integrating sensitivity equations (refs. 7 and 8). The purpose of this example 
is to compare the relative performance of the two approaches on a problem which is 
characteristic of dynamic systems, such as aircraft estimation problems. The impor- 
tant characteristics of the aircraft estimation problem are the requirement to inte- 
grate equations of motion to evaluate the cost function and the requirement to have 
reasonable initial values for the estimated parameters. These initial values might 
be given by a least-squares procedure. 

The first test problem, referred to as estimation problem I, is a single- 
input/double-output, linear second-order system with six unknown parameters. The 
system equations used to generate simulated data are given by equations (1) and (2). 
The six unknown parameters are the four elements of the 2x2 system matrix A and 
the two elements of the control-input matrix B found in equation (1). Matrices A 
and B are defined in table I as true values of 0; matrices C and D of equa- 
tion (2) are the unity and zero matrix, respectively. The data were generated 
by using a simple Euler integration technique, and the input was selected as 
U = sin(t). Initial conditions were set to 0. Process noise and measurement noise 
are excluded for problem I. 

For aircraft estimation problems, the bulk of computer time is spent in perform- 
ing integrations. To prevent any bias in the results due to variations in program- 
ming efficiency or integration techniques, the estimation algorithms use the same 
integration subroutine. For estimation problems I and II a simple Euler integration 
technique is used and the integration subroutine is specifically designed to inte- 
grate matrix second-order linear time-invariant systems. In order to accommodate 
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TABLE I.- ESTIMATION PROBLEM I 


Unknown parameter, 
0 

True 

value of 

Initial 
value of 

Final estimated values 
using method - 

0 

0 

SM 

MNR 

MNRES 
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At, sec 
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the MNR algorithm, which requires a fourteenth-order system to be integrated each 
pass, the system 
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is treated as seven, second-order systems: one second-order system for the states 

(given by eq. (1)), four second-order systems for the sensitivities with respect to 
A matrix elements (given by eq. (13) with j = 1 through 4) and two-second order 
systems for the sensitivities with respect to B matrix elements (given by eq. (14) 
with j = 1 and 2). The results of the first-system integration X(t) are used as 
the input function for the X^ system. Because the systems are not fully coupled, it 
is more efficient and just as accurate to integrate them separately (as described) as 
it is to integrate the complete fourteenth-order system. 

The performance of the methods used in this report will be evaluated with the 
following criteria: 

(1) Accuracy of estimates 

(2) CPU time to termination 

(3) Number of cost-function evaluations 

The primary evaluation criteria will be accuracy of estimates and CPU time to termi- 
nation. Termination is obtained when parameter and cost-function fractional changes 
are computed to be within a specified precision. Both cost-function change AJ/J 
and parameter change A0/9 are required to be satisfied simultaneously to prevent 




premature termination on a plateau where AJ << 1 and A0 is relatively large, or on 
a steep slope where A0 << 1 and AJ is relatively large. 

A secondary criterion is the number of cost-function evaluations, which is used 
only as an indirect guide since each method advances differently with each cost deter- 
mination. However, the number of cost-function evaluations is important since each 
method requires the system equations to be integrated at least once for each cost 
evaluation, with one cost evaluation representing a significant computational effort. 
To emphasize that the same measure of computational effort is used to compare the 
different methods, the term "equivalent evaluations" is used. One equivalent evalua- 
tion represents the computational effort required to integrate the system equations 
once. Each method described in this report requires a different number of equivalent 
evaluations to make one update in the parameter estimates. 

Only dynamic systems or aircraft system estimation problems are used in this 
study, rather than classical test problems such as Rosenbrock's function (ref. 3). 
Using classical optimization problems, which usually require very little computational 
time to evaluate the cost function, could lead to different conclusions about the 
algorithms. The problems in this study range from a simulated linear system with 
6 unknown parameters and noise-free measurements (problem I) to a real aircraft 
problem modeled by a nonlinear system with 20 unknown parameters and measurements 
corrupted by noise (problem IV). 

Figure 1 shows the input and response time histories for problem I, and table I 
gives the numerical results. This table shows the true value of the six unknown 



0 10 20 
t, sec 


Figure 1 .- Time history of input and 
response variables for problem I. 
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parameters, their initial values, and the final estimated values. Each method was 
able to estimate the parameters correctly and satisfy the convergence criterion. The 
MNR method performed about 30 times faster than the SM. This is apparent when con- 
sidering the number of equivalent evaluations required by each method, 715 for SM and 
28 for MNR. 

Table I shows clearly that optimization problems having reasonable starting 
values and involving time-consuming cost-function evaluations should not be solved 
with direct search methods, such as SM. Reasonable initial values tend to provide a 
more quadraticlike cost function for which Newton's method is most effective. If 
reasonable initial values are not available, the SM may be more attractive. In 
response to the results of problem I, this study concentrated on improving the MNR 
method. The next section describes the optimization method developed from MNR. 


MODIFIED NEWTON-RAPHSON METHOD WITH ESTIMATED SENSITIVITIES ( MNRES ) 

The MNRES method developed in this paper is essentially an MNR optimization algo- 
rithm with an efficient method for estimating sensitivities. The sensitivities are 
determined with less computational effort than by using either a finite-difference 
method or analytically determined sensitivity equations. As in the ML/MNR algorithm 
previously described, the same equations (eqs. (1) through (12)) apply for ML /MNRES ; 
however, the sensitivities G. are computed by using slope information from local 
surface approximations of Y(0). The fitted surface allows the slope or sensitivity 
information to be updated at each iteration with a significant reduction in compu- 
tational effort. 

In reference 10, a nonlinear least-squares algorithm is presented which uses a 
surface-fitting method. This algorithm uses a linear-surface approximation of a 
scalar-response variable to eliminate derivative calculations altogether. In MNRES, 
the surface approximation is treated differently by developing an algorithm which 
retains derivative information in a Newton— Raphson method for multivariable systems. 
This is done to provide directional information for the convergence process and to 
provide covariance information. 

The sensitivities for MNRES are determined by assuming that they are approxi- 
mately equal to the slopes of a surface which has been fit to Y(9). The surface 
approximates Y(0) near the series expansion point of equation (8). The choice of 
the type of surface (e.g., nth-order polynomial or spline) and the method of fitting 
the surface (e.g., least squares or simply solving n simultaneous equations for 
n unknowns) is made by the user to suit the particular need. The trade-off in choos- 
ing a method involves the choice between accuracy of the sensitivities and computa- 
tional effort. The simplest and computationally least demanding approach is to 
expand Y(9) in an n-term first-degree polynomial in 0 and then solve for the coef- 
ficients (sensitivities) by solving the n simultaneous equations. Although this 
approach produces the least accurate sensitivities, this, as well as any other 
approach, can have the improved accuracy through step-size control, if necessary. 
Theoretically, in the limit as step size becomes very small, the estimates from the 
surface fits would equal the estimates from a finite-difference method, which in the 
limit equals the actual slope. The trade-off for step-size control involves the speed 
of convergence versus the accuracy of the sensitivities and robustness of the algo- 
rithm to noisy data. 
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The MNRES algorithm is best described by first looking at the computationally 
least demanding approach of using a linear-surface approximation. Expanding Y(9) 
in a first-degree polynomial in 9 for each point in time and at n + 1 different 
points in the n parameter space gives 


4 (S3) " s k0 


+ + 


+ s 0 J 
kn n 


(16) 


where i indicates the ith point in time? k indicates the kth element of the output 
vector Y( 9 ) ; and j indicates one <?f the n + 1 sample points used to fit equa- 
tion (16) to Y(0). Note that' Y^ (0^ ) = Y(0) at each of the n + 1 points. The 
sample points are chosen by allowing a small perturbation of each parameter. Alter- 
natively, the perturbation size can be selected to reflect the relative significance 
of each parameter to the model. This allows for larger perturbations of the less 
sensitive parameters and smaller perturbations for the very sensitive parameters, thus 
providing higher quality derivative calculations. This alternative is discussed 
further in the next section. The slopes s fc1 to s kn are the desired sensitivities 
{Sy^/dQ •} j_f and s kQ is the point on Y(0) where the sensitivities are desired; it 
is the series expansion point of equation (8) for the MNR optimization method. Note 
that because this is a linear surface, the slopes are constant over the surface and 
need not be evaluated specifically at s k Q* a higher degree polynomial is fit to 

Y (0 ) , the slopes will vary across the fitted surface and, therefore, must be evaluated 
specifically at s kQ . Consider the matrix representation of equation (16) for the 
first element of Y and for the n + 1 sample points: 
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Since s 1Q is a known point, equation (18) can be simplified. The first line in 
equation (18) can be eliminated by subtracting from the other n equations. Thus, 


AY, . = AX S 
lx 1i 


( 19 ) 
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( 20 ) 


Thus, at time i the sensitivities for the first element in Y are given by 


S = [AX]” 1 AY,. 
1 1 1 1 


(21 ) 


Note that the AX matrix is independent of time. This enables the sensitivities 
to be calculated rapidly during each iteration of the algorithm. This is a key 
factor in reducing the computational effort of the algorithm; in effect, the integra- 
tion of the mn sensitivity equations has been replaced by a set of m matrix 
multiplications. 

Figure 2 shows, geometrically, two iterations for the case in which 0 is of 
dimension two and a linear surface is used to fit a scalar y. The expansion is 
simply at time i: 



i 


( 22 ) 


During the first iteration, this espansion requires that y(0) be evaluated at 
n + 1 = 3 points: y° , y 1 , y 2 . Computationally, the first iteration is the most 

costly phase of the MNRES algorithm. Each evaluation of Y requires that the equa- 
tions of motion be integrated. The linear surface (indicated by the solid-line tri- 
angle in fig. 2) is fit and the slopes (sensitivities) are thereby determined. The 
algorithm proceeds as in the ordinary MNR method to obtain 
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Figure 2.- Linear-surface fit for two iterations of MNRES. 


3 

The new Y is evaluated (through integration) at 0 k+1 to 9 et Y O)* At this point 
the MNRES algorithm has reduced the sensitivity problem to solving a set of simulta- 
neous equations. This is done by eliminating the 9-^ in X, which produced the 
greatest cost in J(0), and replacing that information with the newest estimate 
of 9. The new surface (indicated by the dotted-line triangle) in figure 2 assumes 
that y° is the high-cost point and thus eliminates it from the fitted surface. The 
slopes of the new surface provide the sensitivities for the MNRES algorithm to pro- 
ceed. In this method a check should be made to ensure that the new Y-^(9 j, + ^ ) pro- 
duces a smaller value of J(0). In some cases, step-size control or complete 
restarting may be needed. Note that initialization of the algorithm requires that 
n + 1 integrations be performed for the n + 1 trajectories, Y^ • After the 
start-up, only one integration of the system equations is needed to evaluate the 
cost J(9), outputs Y( 0 ) , and updated parameter estimates for each iteration. 

The least-squares approach to fitting the surface Y(0) offers another advantage 
if an iterative least-squares method is used. The iterative method provides a memory 
device reducing the storage requirements from n + 1 sets of output time histories to 
just two time histories. One of the two corresponds to the new response predicted by 
the most recent estimate of 0, and the other corresponds to the outgoing 0 that 
produced the highest cost. The penalty for this advantage is the need to integrate 
equations of motion twice per iteration; this result still requires substantially less 
computational effort than that required with the usual MNR method. This algorithm was 
used in estimation problems I and II and is described subsequently and also in the 
appendix. It should be noted that because problems I and II did not require a large 
amount of storage space, all time histories were saved so that only one integration 
per pass was performed . 

When using the iterative least-squares approach, only two changes are made to 
the MNRES algorithm just described. The first change is in the calculation of the 
AX matrix, and the second change is in the sensitivity calculation. The development 
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of this formulation will begin with equation (16) in condensed form, which is shown in 
the following equation. Everything discussed up to this equation in the previous 
development applies here. Thus, 



X : S. 


ki 


(23) 


Simplifying the notation by dropping the ki subscript and writing the matrix form of 
the equation (which removes the j superscript) gives 


Y = XS 


(24) 


The least-squares solution for the sensitivity vector is 


S 


T -1 T 

[XX] 'x Y 


(25) 


Now, defining a recursive relation for the k + 1 iteration gives 


Vi - K )T \T 


(26) 


and the updating equation is 


' k+1 


-1 . , .T / o\T O 

p k + la k> \ - (“kj “k 


-1 


(27) 


where is the new set of 0 to be included and is the outgoing set of 0 

’ “ "* ... )c , . 


which produced the high cost in J. The recursive equation for S is 


s k*i - \ - \ +1 


/ o\T O , .T 

(\) y k - l \> y k 


T 

( V \ 


- 


(28) 


With the new sensitivities determined, the algorithm proceeds as before. The recur- 
sive relations are derived in the appendix. 


PROPERTIES OF MNRES OPTIMIZATION METHOD 

In this section, properties of MNRES which aid convergence are discussed heuris— 
an< ^ factors which have a significant effect on speed are indicated. Also, 
similarities and differences with the commonly used MNR method are discussed providing 
a benchmark for comparison. Finally, to demonstrate the performance of MNRES, two 
estimation problems with varying noise levels are solved by using MNRES and MNR to 
obtain ML parameter estimates. 
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The ML parameter estimates are obtained by solving an unconstrained, nonlinear 
optimization problem; that is, find 9* which minimizes the cost function J(9). The 
necessary and sufficient conditions for this problem are as follows: 

(1) J( 0) is differentiable at 0*. 

(2) VJ(0*) = 0. 

(3) V 2 J(0*) > 0. 

The theory for solving the unconstrained optimization problem is often based on the 
assumption that the cost function J(0) is a quadratic function of 0. This provides 
a much more tractable theory and allows basic properties to be readily established. 
Corresponding theorems for general nonlinear functions are very difficult to prove. 
However, techniques developed using the quadratic assumption are still very effective 
for nonlinear functions. Many techniques for solving nonlinear minimization problems 
are developed from practical experience. 


Comparison of MNR and MNRES Properties 

Convergence of NR or MNR algorithms, both with and without finite-difference 
derivatives, has been well documented (ref. 3) • Convergence of MNRES can be shown, at 
least heuristically, by considering several details. First, the MNRES method is still 
fundamentally an NR method or, for this study, an MNR method. The only critical dif- 
ference is that the derivatives are approximate which makes MNRES closer to MNR with 
numerically determined derivatives. Second, note that fitting a first-degree, n-term 
polynomial such as 


y 


s + s 0 + 

0 1 1 


+ s 0 
n n 


( 29 ) 


to n + 1 data points is equivalent to a simple finite-difference method. In effect, 
as A0^ (the distance between points on the fitted surface for MNRES) becomes 
small enough, the sensitivities become identical to that given by a simple finite- 
difference method, regardless of the actual functional representation of Y(0). The 
MNRES algorithm simply relaxes the accuracy of the sensitivities in order to reduce 
substantially the integration requirements; the degree of relaxation varies during the 
optimization process but can be controlled by limiting step size. 

The relaxation of sensitivity accuracy generally appears to be a very beneficial 
trade-off for Newton-Raphson algorithms . During an MNRES optimization there are two 
times that the sensitivities are very close to the finite-difference values; these 
times are during the initialization or start— up of the algorithm and toward the end 
when A9 becomes small. During initialization of the algorithm, n different 0 
are chosen (a perturbation on each element of 0 is sufficient) and the actual sur- 
face given by Y(0) is fitted. The user can choose an initial 03 such that 

1(03 _ 0 ° ) /e° | « 1 for each j. In this study, the same A0^ was used in the MNR 
with finite-difference derivatives as that used in the start-up of the MNRES. This 
was done for comparison purposes; in practice, the choice of perturbation size for 
0 may be very different, as discussed later. 
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The computational advantage of MNRES over MNR is determined to a large extent by 
the number of unknown parameters n. During the first iteration, both MNRES and MNR 
require n + 1 equivalent evaluations. After the first iteration, MNR continues to 
use n + 1 equivalent evaluations per iteration but MNRES needs only one per itera- 
tion. As the estimation process advances, MNRES continually eliminates values of 0 
which are far from 0*, based upon the information in one equivalent .evaluation. 

Thus, at the end, 0 is very close to 0*, producing very small A0-* for the deriv- 
ative estimations. Again, the derivatives obtained by the surface-fitting method 
approach the values obtained in a finite-difference method. 

Two factors aid in preventing divergence during the critical time period, when 
MNRES is between initialization and the end of the optimization, in which potentially 
large A0 can occur. First, the user can always incorporate step-size control logic. 
Carried to the extreme, MNRES could always be forced to approximate the derivatives 
the same as a finite-difference method. Of course, convergence would be very slow 
because of the very small steps. In practice, one would let the algorithm take steps 
determined by the NR logic (as done in this study); and then if a convergence problem 
develops, one would begin controlling step size. Secondly, NR, MNR, and MNRES methods 
advance more quickly as the quadratic approximation of the cost function improves; 
moreover, the Newton algorithm converges in one step for a quadratic cost function. 
Since the quadratic approximation of the cost function improves the closer that 0 
gets to 0*, and since initial estimates of 0* are often given by a least-squares 
procedure or knowledgeable user, 0 q tends to be close to 0*. Thus, MNRES starts in 
a region conducive for convergence and, therefore, the method would not be required to 
handle very large A0 for many practical cases. 


As mentioned previously, in practice it is beneficial to choose the perturbation 
size in a different fashion from that used in a simple finite-difference method. 

Simply using a 1 -percent perturbation on each element of 0 to obtain the correspond- 
ing perturbation in each element of Y(0) is not optimum for derivative calcula- 
tions. Experience has shown that it is beneficial to use perturbation sizes which 
reflect the importance of the parameter to the model. By computing the sensitivities 


as 


°j"jj 


for each parameter and then letting the perturbation sizes be scaled 


inversely proportional to the normalized ratios of sensitivities, more accurate deriv 
ative information can be obtained. Of course, this applies only when an initializa- 
tion or "restart" is needed. The fundamental issue is that the less sensitive a 
parameter, the larger the perturbation necessary to obtain an appropriate size 
response in the outputs. This approach could also be applied to an MNR method. 
Theoretically, the same derivative should be obtained for any sufficiently small 
perturbation in 0; however, because of both the sometimes widely varying sensi- 
tivities of the parameters and the numerical-precision limitations, it is beneficial 
to vary the perturbation size according to the aforementioned rule. The sensitivity 


defined as ® j M j j was introduced in reference 11 and used again in reference 12 as a 
means of quantifying the significance of a parameter to the model. 


There are four factors which have a significant effect on the speed of the MNRES 
algorithm. The first factor is the degree of precision demanded by the convergence 
criteria AJ/J and A9/0. Naturally, the more precision desired, the greater the 
computational demand. How well the algorithm meets the demand is determined by the 
computer capability (word size) and by the information content of the data analyzed. 
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The second factor is the accuracy of the sensitivities, particularly during initial- 
ization of the algorithm. The perturbation size chosen for 0 can have a significant 
effect on the accuracy of the sensitivities. In turn, the accuracy of the sensitivi- 
ties directly control the speed of the algorithm since G. determines both the step 

size and direction in the parameter-update equation. A third important factor is the 
scaling term used in most ML algorithms to reduce the step size. This term is com- 
monly used to aid convergence since the computed step size is sometimes too large and 

may cause divergence. The scaling term used in this study was equal to 1 • The fourth 
factor is unique to MNRES and is determined by the number of parameters n to be 
estimated. This problem is a result of inverting the n x n matrix AX, which may 
contain very small numbers as the algorithm proceeds. Experience has shown that 
checking this matrix inversion for numerical difficulties is important in the MNRES 
method. 


Problem I Extended 

Estimation problem I is solved again to allow a comparison between MNRES and the 
two methods previously considered, MNR and SM. Problem I has been described and fig- 
ure 1 shows the input and response time histories. MNRES uses the same finite- 
difference method as that of MNR to determine sensitivities during initialization. 
Also, MNRES uses the iterative least-squares form of the algorithm; however, only one 
integration per pass is performed. Because of the small storage requirements, all 
time histories are saved, thus eliminating the need for an additional integration. 
Table I presents the results of using the three optimization methods in the estimation 
algorithm. All three algorithms accurately converged to the correct parameter values. 
However, MNR was 30 times faster than SM; and MNRES was twice as fast as MNR, or 
60 times faster than SM. The number of equivalent evaluations had similar ratios, 
that is , 71 5 : 28 : 1 2 . 

The initial values were relatively close to the final solution and thus allowed a 
good quadratic approximation of the cost function; therefore, the two Newton-Raphson 
methods were substantially faster. MNRES, however, capitalized more efficiently on 
the information obtained from each integration of the system equations. Each integra- 
tion of the system equations provides information which is immediately incorporated 
into the numerical process when using MNRES. When using MNR, n + 1 system integra- 
tions (equivalent evaluations) are required before each updating operation; for prob- 
lem I, n + 1 is equal to 7. 


Performance of Methods on Estimation Problem II 

Estimation problem II is provided to demonstrate the robustness characteristics 
of MNRES compared with those of the commonly used MNR. The form of MNR and MNRES used 
in problem I is used again in problem II. The known system from problem I is analyzed 
again in problem II, except that measurement noise is added and a step input is used 
to excite the system. Two cases are considered with different levels of measurement 
noise. The noise is zero mean and Gaussian; the standard error of the noise for each 
case is 0.0001 and 0.001, respectively. Figure 3 shows time histories of the input 
and response variables for the two cases, and table II shows the estimation results. 

In case 1 both methods produce equally degraded results; however, MNRES still 
converges to the same precision level more quickly. In case 2, by using a severe 
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noise level and a limited amount of data, MNRES was unable to converge at all. The 
results showed that it was oscillating about a solution, unable to find a new param- 
eter vector which would produce a lower cost. The MNRES used on this problem had no 
special step-size control logic. The solution that was obtained, however, was as 
accurate as that obtained by MNR, which did converge. 
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(a) Case 1* 

3.- Time history of 


(b) Case 2. 

input and response variables for problem II • 


Meeting convergence requirements does not guarantee accurate results; the error 
in the estimates ranged from a 5-percent error to a 130-percent error. MNR had both 
the most accurate and the least accurate estimate. The importance (sensitivity) of a 
parameter to the model will significantly affect the accuracy of the estimate, partic- 
ularly under these adverse conditions. Based on these examples, it appears that MNRES 
performs faster than MNR while providing the same level of precision. 

It should be noted that a Euler integration method was used and may have a 

slightly degraded performance, particularly for the MNR method. Euler integrations 

tend to have an increasing error with the length of integration. This problem was 
reduced in part by using relatively few data points. In problem I, only the first 
5 sec of data (20 data points) were analyzed. In problem II, for both cases, only 

3 sec of data (12 data points) were used. Additionally, step-size control of the 

optimization algorithm could have been used as another means of improving the con- 
vergence capability. This was not used by any algorithm, although SM has the ability 
to determine whether to extend or contract the polyhedron. 
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TABLE II.- ESTIMATION PROBLEM II 


(a) Case 1 


Unknown parameter, 
0 

True 

value of 

Initial 
value of 

Final estimated values 
using method - 

9 

e 

MNR 

MNRES 

CD CD 

o — > 

0 

0.01 

-0.0675 

-0.0684 

-1.5 

-1.6 

-1.471 

-1.471 

6 3 

1.0 

i.i 

1.009 

1.010 

0 4 

-.5 

-.6 

-.449 

-.448 

% 

.2 

.25 

.202 

.202 

% 

.1 

.15 

.098 

.098 

Cost 

Equivalent evaluation • • • 
At, sec 



0.105 E-06 
42 
77.54 

0.105 E-06 
1 2 
24.08 


(b) Case 2 


Unknown parameter, 

True 

value of 

Initial 
value of 

Final estimated values 
using method - 

0 

0 

9 

MNR 

1 

MNRES* 

9 1 

®2 

0 

0.01 

-0.705 

-0.410 

-1.5 

-1.6 

-1.228 

-1.549 

0 3 

1.0 

1.1 

1.757 

.799 

®4 

-.5 

-.6 

.1 59 

-.238 

0 5 

.2 

.25 

.210 

.251 

06 

.1 

.15 

.087 

.037 

Cost 

Equivalent evaluation • • • 
At, sec 



0.104 E-04 
70 

134.44 

0.122 E-04 
27 
56.05 


^ MNRES did not converge. 


APPLICATION TO NONLINEAR AIRCRAFT ESTIMATION PROBLEMS 

In this section the utility of MNRES is demonstrated by application to nonlinear 
aircraft estimation problems by using both simulated- and real— flight data. For 
these problems, the computationally least demanding form of MNRES (with eq. (21) 
instead of (28)) is used to compute sensitivities. The accuracy and robustness of 
the algorithm are assessed by testing the algorithm on simulated data with known 
noise levels. Also, performance of the algorithm on real-flight data is checked by 
comparison with a well-proven ML algorithm. In both cases the problems involve 
nonlinear lateral models of a general aviation aircraft. 
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Estimation Using Simulated-Flight Data 


In this section ML estimation problem III is considered. This problem simulates 
real-flight data with varying noise levels. Three cases are considered: case 1 

without any measurement noise, case 2 with a representative noise level typical of 
flight data for the aircraft, and case 3 with twice the noise level of case 2. The 
standard errors of the simulated measurement noise are shown in table III. In each 
case the noise is zero mean and Gaussian. The simulated data were created by a 

TABLE III.- STANDARD ERRORS OF SIMULATED 
MEASUREMENT NOISE 


Output variable 

Standard deviations for - 

Case 1 

Case 2 

Case 3 

[3, rad 

0 

0.01 0 

0.02 

p, rad/sec .... 

0 

0.010 

0.02 

r, rad/sec .... 

0 

0.01 0 

0.02 

4 ), rad 

0 

0.005 

0.01 

a^, g units ... 

0 

0.005 

0.01 


fourth-order Runge-Kutta integration with a step size of 0.05 sec. A maximum 
likelihood program called MAX is used to solve the problem and it uses the MNRES 
optimization method. The MNRES method computes sensitivities by using a simple 
finite-difference technique during initialization and then uses equation (21 ) there- 
after. Program MAX uses two convergence criteria: AJ/J < 0.001 and A0/0 < 0.001 • 

Table IV shows the terms used in the nonlinear aerodynamic model to create the 
simulation and the parameter estimates obtained through analysis of the simulated 
data. Time histories are provided for the three cases in figure 4. The control 
inputs were the same for all three cases and are shown in figure 5. As expected, the 
estimates of the less easily identified terms are more quickly corrupted as the noise 
levels increase; however, the estimates are still very reasonable and the time histo- 
ries are accurately predicted. Table IV shows that the MNRES method can be used 
effectively in estimating parameters for nonlinear aircraft systems. 


Estimation Using Real-Flight Data 

In this section ML estimation problem IV is considered. Problem IV uses flight 
data from a general aviation aircraft operating at an angle of attack of 8°. The 
problem is solved by using two different ML programs. The first, referred to as 
program MAX, uses the MNRES algorithm; and the second, referred to as program MAXLIK, 
uses an MNR algorithm. This MNR algorithm integrates the sensitivity equations to 
obtain the sensitivities. As used in the simulated-data case, the MNRES method com- 
putes sensitivities by using a simple finite-difference technique during initializa- 
tion and then uses equation (21) thereafter. MAXLIK is a proven code used at the 
Langley Research Center for aircraft-parameter estimation. MAXLIK was developed from 
the equations documented in reference 5. For comparison purposes each program uses a 
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TABLE IV.- SIMULATED-DATA ANALYSIS USING MNRES 


Unknown 


Parameter 

estimates 

for - 

parame tar r 

Simulation values 




e 


Case 1 

Case 2 

Case 3 

O 

>* 

u 

0.1 3 

0.1 299 

0.1 298 

0.1 295 


-.411 

-.4136 

-.4261 

-.4401 

% 

-.146 

-.1 524 

-.1 874 

-.2379 

c Y 

r 

.63 

• 6686 

.6070 

.5412 

C Y 

Y 6a 

-.053 

-.0618 

-.0733 

-.0872 

c v 

Y 6r 

.075 

.0794 

.0775 

.0751 

C i.O 

0 

.0001 

-.0003 

-.0005 


-.123 

-.1223 

-.1228 

-.1 240 

c x 

-.397 

-.3988 

-.4026 

-.4094 






C x 

.257 

.2573 

.2409 

.2239 

x r 





C x 6a 

-.1 82 

-.1 815 

-.1778 

-.1 755 

° l 6t 

.077 

.0067 

.0059 

.00497 

C. 

2.63 

2.6254 

2.519 

2.4359 

op 





O 

3 

O 

0 

-.00005 

-.00008 

-.0001 

° n P 

0 

.000003 

.0001 

.0005 

% 

-.15 

-.1488 

-.1 524 

-.1558 

Cn r 

-.083 

-.0828 

-.0861 

-.0911 

Cn 6a 

-.0431 

-.0425 

-.0434 

-.0445 

Cn 6r 

1.7 

1.7343 

1.4419 

1.0118 
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Figure 5.- Control inputs for 
lateral simulation. 


fourth-order Runge-Kutta integration method and an integration step size of 
0*05 sec* A convergence criterion is set at AJ/J = 0*001 for both codes. Pro- 
gram MAX normally uses an additional criterion restricting the parameter change 
A9/0; however, in this problem it is removed to ensure that both programs converge 
for the same criterion. Both programs use the same bias and scale-factor corrections 
to the flight data. These corrections were determined by using a compatibility 
program developed in reference 12. 

Estimation problem IV involves a nonlinear lateral model. Table V presents a 
comparison between parameter estimates and their standard errors from both MAX and 
MAXLIK. It also shows the initial values from the modified stepwise regression 

2 

program of references 1 and 13. The sensitivities computed as 0 j M jj are 9 lven in 

the last column. Again, there is reasonable agreement between the two approaches. 
Standard-error estimates tend to be a little higher for program MAX; this is probably 
due to their sensitivity to the derivative information. Repeating the calculations 
with program MAX, by allowing the sensitivity ratios to be incorporated into the 
initializing derivative calculations, provided a small improvement in the overall 
speed of the algorithm. This occurred because only one restart was required during 
the optimization process. More improvement would be realized in problems where 
restarting occurs several times. Time histories of the measured flight data and 
predicted response using the estimated model are shown in figure 6. Execution times 
for problem IV indicate program MAX to be about 30 percent faster. 
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TABLE V. - ESTIMATION PROBLEM IV 


Unknown 

parameter, 

9 

Initial value 

Program MAX 

Program 

MAXLIK 

0 2 M.. 
D 33 

of e 

a 

0 

A 

a 

A 

6 

A 

a 

c 

Y,0 

0.036 

0.0061 

0.0006 

0.0213 

0.0005 

0.7533 E+03 

% 

-.479 

-.4603 

.0075 

-.4608 

.0067 

.3756 E+05 

° Y P 

-.186 

-.1378 

.0485 

-.0604 

.0439 

.9373 E+03 

C Y 

r 

.522 

.6677 

.0289 

.6209 

.0256 

.1915 E+04 

C Y 

Y 6a 

-.08 

-.0504 

.0166 

-.0375 

.0150 

.8129 E+03 

c v 

Y 6r 

.083 

.0814 

.0043 

.0763 

.0037 

.9922 E+03 

c Y 

Y ap 

.45 

.4300 

.0592 

.4512 

.0504 

.9399 E+03 

°x,0 

0 

.0002 

.00005 

-.0001 

.00005 

.9618 E+03 

° l P 

-.079 

-.0872 

.0015 

-.08 

.0013 

.1493 E+06 

C ' P 

-.47 

-.5320 

.0102 

-.4823 

.0085 

.4529 E+07 

% 

.187 

.1700 

.0043 

.1543 

.0045 

.6838 E+04 

Sa 

-.19 

-.2035 

.0036 

-.1852 

.0031 

.7969 E+05 

° x 6r 

.01 

.00055 

.00024 

-.0012 

.00072 

.1340 E+04 

C 'ctp 

-.26 

-.2707 

.0116 

-.2105 

.0091 

.6327 E+04 

C n, 0 

0 

-.00063 

.00003 

-.002 

.00002 

.1909 E+04 


.04 

.0323 

.00045 

.0329 

.0004 

.1515 E+06 

% 

-.056 

-.1043 

.0026 

-.0916 

.0022 

.1400 E+06 

C "r 

-.15 

-.1462 

.002 

-.1534 

.0017 

.5039 E+05 


’ 0 

-.0044 

.001 

-.0037 

j 

.0009 

.1780 E+04 

o 

3 

o> 

-.053 

-.0550 

.0003 

-.0532 

.0003 

.9048 E+07 

V 

V 

-.39 

-.39 


-.39 



1c 3 
n p 3 

.08 

.08 


.08 




Parameter held fixed. 
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Figure 6.- Real-flight data from measured and predicted 
responses for lateral simulation. 
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Figure 6.- Concluded. 


CONCLUDING REMARKS 

An algorithm for maximum likelihood (ML) estimation using an efficient method 
for estimating the sensitivities has been developed. The optimization method is 
referred to as a modified Newton-Raphson with estimated sensitivities (MNRES) • The 
method is applicable to any parameter-estimation problem and is particularly suited 
for multivariable dynamic systems. In this study nonlinear aircraft systems were 
analyzed. 

The algorithm has three advantages over other commonly used techniques. The 
first advantage is that the algorithm removes the need to derive sensitivity equa- 
tions for each new model; this eliminates the computational burden of integrating the 
sensitivity equations during each iteration of the algorithm. This also provides a 
lot of flexibility, allowing the model equations to be in any format that is con- 
venient - such as splines, polynomials, or a nonanalytic form. The second advantage 
is that the algorithm is effective for a variety of methods chosen to fit the output 
vector surface in parameter space (needed for sensitivity estimation), allowing the 
user to choose a surface-fitting method best suited to the problem. An approach is 
discussed which provides a method to reduce storage requirements with little addi- 
tional computation. The third advantage of the algorithm is that it reduces the 
computational effort in comparison with the commonly used approach referred to as a 
modified Newton-Raphson (MNR) method. For small problems (fewer than 15 parameters 
to be estimated), the reduction can be substantial. For larger nonlinear problems, 
the reduction is usually more modest; however, improvements may still be significant 
if data quality, signal compatibility, and sensitivity calculations are good. Based 
on this study, it appears that the ML /MNRES algorithm generally performs better than 
the commonly used ML /MNR algorithm. 

Langley Research Center 

National Aeronautics and Space Administration 
Hampton, VA 23665 
May 11, 1984 
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APPENDIX 


ITERATIVE LEAST-SQUARES FORMULATION FOR MNRES 
Consider the least-squares problem 
Y = XS 

and the least-squares solution 

- r T *.-1 T 
S = [x xj X Y 

Define for the i + 1 iteration 



(A1 ) 


(A2) 


(A3) 


Let a_^ be one row of X containing old information to be removed from X, and let 
a^ be a replacement row containing new information to be added to X. Define Z as 
the common elements of X between two iterations. Partitioning X for the i - 1 
and ith iterations results in 


x i-i 


(A4) 


X. = 
1 


a. 

i 


Z. 

1 


(A5) 


By using equations (A4) and (A5), the following relations can be written: 


<vA-i ■ * <v\ 


(A6) 


T T T 

(X i ) X jL = (a^ &i + (Z i ) Z ± 


(A7) 


From equation (A6) obtain 


<V\ - «*!.,» V, - (-?)’-? 


(A8) 
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Substituting equation (A8) into (A7) gives 


<v T *i - - K) T *° - 


Substituting equation (A9) into (A3) gives 


- [ (X i-, ,T *i-, ' K) T »i * <«!>%]"’ 


which can also be written as 


[ -1 r o\T o V T ”1 

p t - K) a i + (a i> a iJ 


Applying the same development to equation (A2) gives 


S . = P. 

i+1 i+1 




and substituting equations (A2) and (A3) delayed a step into (A12) 


■ p i +1 [ p l' s i - ( a I) Tp ° + (a i )T yi] 


Expanding equation (A13) gives 


s i+i = p i + i p i s i - p i+ i( a °) y°i + p i+ i (a i )T yi 


Noting that 


S. - P. p.'s. = 0 
l 1+1 1+1 i 


and then adding equation (A15) to (A14) gives 


s i*i - s i - p i« p I + , s i + p 1M p I - p it ,(»i) T ^ * p i+ ,< a i> T yi 

s i« - s i - p i + ,[ p u, 9 i - p l\ + ( a t) Tp i - <a i )T yJ 


(A9) 

(A10) 
(All ) 
(A12) 
(A13) 

(A14) 

(A15) 

(A16) 

(A17) 
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Combining terms gives 


s i + i - s 


i - p i + i[( p iii - p i i ) s i + ( a i) Ty ° - u i> Ty i] 


and using equation (All) yields the desired relation 


S. = S. 

1+1 l 


i - p i + i[[ <a t )Ta i - ( a i) Ta i] s i + K) Ty i - <a i ,Ty i] 
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